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Abstract 



We have studied the superfluid density p s on various size-lattices in the geometry L x L x H by 
numerical simulation of the x — y model using the Cluster Monte Carlo method. Applying the Kosterlitz- 
Thouless-Nelson renormalization group equations for the superfluid density we have been able to extrap- 
olate to the L — > oo limit for a given value of H. In the superfluid phase we find that the superfluid 
density faithfully obeys the expected scaling law with H, using the experimental value for the critical 
exponent v = 0.6705. For the sizes of film thickness studied here the critical temperature T c and the 
coefficient b entering the equation T/(p s H) cx I — b(l — T/Tc) 1 ^ 2 are in agreement with the expected 
Independence deduced from general scaling ideas. 

1 Introduction 

Liquid 4 He is an ideal experimental testing ground for the theory of phase transitions and the related finite- 
size scaling (FSS) theory. Relevant physical quantities such as the specific heat c or the superfluid density 
p s can be measured to a very high accuracy [l]-[5] and they can be used to check the FSS theory (see e.g. 
[6]). This theory was developed in order to account for the influence of the finite extent of systems which are 
confined in a finite geometry (e.g. a film) and at temperatures close to the critical temperature. The theory 
is based on the scaling hypothesis and on that finite-size effects can be observed when the bulk correlation 
length £ becomes of the order of the relevant size of the system (e.g., in a film the relevant size is the film 
thickness H). More precisely, the finite-size scaling hypothesis states that a dimensionless physical quantity 
(or the ratio of two physical quantities of the same dimensions), sufficiently close to the critical point, is 
a function only of the ratio H/£. For a physical quantity O this simple but non-trivial statement can be 
expressed as follows [7]: 



t is the reduced temperature and / is a universal function. So far the validity of this approach has been 
confirmed by experiments on superfluid helium on Helium films of finite thickness [2] (the relevant size is 
the thickness). However, recent measurements of the superfluid density by Rhee, Gasparini, and Bishop [3] 
of Helium films seem to be in contradiction to the FSS-theory. 

The singular behavior in the thermodynamic functions of liquid 4 He close to the superfluid transition 
can be understood in terms of a complex order parameter ip(r) which is the ensemble average of the helium 
atom boson creation operator. This ensemble average is defined inside a volume of a size much greater than 
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the interatomic distance but much smaller than the temperature-dependent coherence length. In order to 
describe the physics at longer length scales, which is important very close to the critical point, we need 
to consider spatial fluctuations of the order parameter. These fluctuations can be taken into account by 
assigning a Landau-Ginzburg free energy functional 7i(ip(f)) to each configuration of t/>(r) and performing 
the sum of e ~' H / k BT over sucn configurations. The power laws governing the long distance behavior of the 
correlation functions and the critical exponents associated with the singular behavior of the thermodynamic 
quantities close to the critical point are insensitive to the precise functional form of H[ip], and they are the 
same for an entire class of such functionals. The planar x — y model belongs to the class of such Landau- 
Ginzburg free-energy functionals [8], and thus can be used to describe the fluctuations of the complex order 
parameter. In the pseudospin notation the x — y model is expressed as 

n = -j^2s i -8 j , (2) 

(i,3) 

where the summation is over all nearest neighbors, s = (cos 9, sin 9) is a two-component vector which is 
constrained to be on the unit circle. The angle 9 corresponds to the phase of the order parameter t/>(r). 

In this paper we investigate the x — y model in a film geometry, i.e. planar dimensions L with L — ► oo 
and a finite thickness H. This study will allow us to examine directly the validity of the FSS theory. 

There has been analytical and numerical work on the pure 3D x — y model. Results of high-temperature- 
series studies can be found in [9], Monte-Carlo simulations were reported in [10]-[13], a renormalization group 
approach based on vortex lines was reported in [14]. In [15] the anisotropic 3D x — y model (J x = J y ^ J z ) 
was studied. A crossover from 3D to 2D behavior was found with respect to the ratio J z /J x . The Villain 
model, which is in the same universality class as the x — y model, was studied in a film geometry in [16] 
where the correlation length in the disordered phase was used to extract the thickness-dependent critical 
temperature. 

In this paper we study the superfluid density or helicity modulus of the x — y model in a film geometry, 
i.e. on an L 2 x H lattice. In a film geometry this model exhibits a crossover from 3D to 2D behavior. In 
the temperature range where the model behaves effectively two-dimensionally we are able to compute the 
values for the helicity modulus in the L — ► oo limit using the Kosterlitz-Thouless-Nelson renormalization 
group equations. This enables us to eliminate the Independence completely and thus to check scaling of the 
helicity modulus with respect to the film thickness H . We also test further consequences of the FSS theory 
which are described in the next section. In section 3 we show how the helicity modulus is computed in our 
model and briefly describe the Monte-Carlo method. Section 4 discusses the Independence of the helicity 
modulus in the temperature range where the model exhibits two-dimensional behavior. In that section we 
also describe the extrapolation procedure to the L = oo limit. In section 5 we study the scaling of the helicity 
modulus with respect to U and in section 6 we investigate the U-dependence of the critical temperature. 
The last section summarizes our results. 



2 The x — y model and finite-size scaling 

The 3D x — y model shows an order-disorder phase transition. Above the bulk critical temperature T\, the 
spin-spin correlation function decays exponentially. The correlation length characterizing this decay grows 
with the temperature T according to 

i(T) <x (^p") " • (3) 

In the ordered phase the correlation function decays according to a power law. There, the role of a relevant 
finite length-scale is played by the transverse correlation length[17], £t ■ £t, up to a constant factor, is 
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proportial to the ratio T/T, where T is the helicity modulus. Fisher's scaling hypothesis implies that 
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Figure 1: The behavior of the 3D x ■ 
temperature. 



y model in a film geometry for a fixed thickness H with respect to the 



If the model is considered in a film geometry, i.e. infinite planar dimensions and a finite thickness H , 
interesting crossover phenomena take place. In Fig.l we show an intuitive picture of the behavior of our 
model with respect to the temperature for a fixed thickness H . Let us start at a temperature far below T\ 
where £t is much smaller than the thickness H . At these temperatures the model behaves as a true 3D 
system. If we raise the temperature in order to approach closer to T\ from below, the correlation length £t 
grows according to (4) until we reach the crossover temperature T~ < T\ where it becomes comparable to 
H . Above this temperaure the behavior of the system crosses over from the 3D to two-dimensional (2D) 
behavior. A further increase in temperature makes the 2D behavior of the system more and more pronounced. 
The system starts "feeling" a 2D critical temperature T^ D . Very close to T^ D , i.e. in the purely 2D regime, 
we have to apply the Kosterlitz-Thouless-Nelson theory in order to explain the behavior of the system. This 
implies that in this regime the correlation length £t does not depend on T according to Eq. (4). Instead, 
the dimensionless ratio K = T/(TH) satisfies the following renormalization group equations [18]: 



dK(l, T) 



4tt 3 j/ 2 (/,T), 



= (2- 7 rA- 1 (/,T))j/(/,T). 



(5) 
(6) 



In y is the chemical potential to create a single vortex, e l denotes the size of the core radius of a vortex. In 



the limit 



all vortices have been integrated out, and T — ► T ( 
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finds [18] 
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where & is a constant. K(l = oo,T) is infinite above T c 
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K(l = oo) exhibits a universal jump at 



T^ D . Above T^ D but still in the purely 2D region, the correlation length grows with T according to the 
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Kosterlitz-Thouless theory, i.e. 



i T (T)Ke^B/(T-T™yl\ 



(8) 



A further increase in temperature results in reaching another crossover temperature T+ , where the correlation 
length £t is comparable to the thickness H . At higher temperatures the model exhibits pure 3D behavior 
again. From this intuitive picture we deduce the following inequality for the crossover temperatures, T r 



2D 



and T: 
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T~ <T? D < T\ <T+. 



(9) 



Within the FSS theory the behavior of our model in the crossover region at temperatures T < T~ can be 
described by a universal scaling function $ which depends only on the ratio H/^t(T), provided H is large 
enough. According to Ambegaokar et al. [19] we have 



K(T, H) 



T 



T(T, H)H 



(10) 



where t = (T\ — T)/T\. In the argument of $ we have replaced £t(T) by its bulk scaling expression (4). 

Interesting conclusions can be drawn if we extend the validity of Eq. (10) up to T 2D [19]. Since 
K(T, H) drops discontinuously to zero at T^ D , the scaling function $ has to be discontinuous too. From the 
Kosterlitz-Thouless-Nelson theory we have to require 
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x c is a dimensionless number. From here we immediately derive an expression for the if-dependence of T c 
namely [19] 
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T? D (H)=T X (l 



(12) 



Eqs. (7) and (10) have to be reconciled in the two-dimensional region. This enables us to deduce the form 
of the universal function $(«) for values of x close to x c . In order to do that we replace T^ D by T\ in Eq. 
(7) by inverting Eq. (12). Keeping only terms linear in H -1 !" we obtain 



K(T, H) 



H l/2v 



1/2- 



For T close to T^ D (i.e. t close to 1 1 — T£ D /T x \) we have 

X r- 



t 



which yields 



K(T, H) 



H l/2v 



tH 1/v 



1/2- 



K(T, H) can only be a universal function of tH 1 !" if b scales according to 

b(H) = AH 1 ' 2 ", 



(13) 
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i.e. we obtain 

= l(l-A(x-x c ) 1/2 ), (17) 

A is a constant which will be found numerically (see later in the text). The form (17) of the universal 
function is valid at temperatures very close to T 2D , since Eq. (7) is only an approximation itself. Petschek 
[20] uses a different argument in order to derive (16). 

3 The helicity modulus and Monte-Carlo method 

The helicity modulus T was introduced by Fisher, Barber and Jasnow [21]. It is related to the superfiuid 
density via [21] 

p s (T) = (jf?(T). (18) 
For the 3D x — y model on a cubic lattice the definition of the helicity modulus is [10, 22]: 

77 \J2 cos ( di - 0j){e^ -Ci/) 2 



T 

1 m 





(19) 

where V is the volume of the lattice, is the unit vector in the corresponding bond direction, and is 
the vector connecting the lattice sites i and j. In the following we will omit the vector index since we will 
always refer to the ^-component of the helicity modulus. Note that, because of isotropy, we have = T y . 
The above thermal averages denoted by the angular brackets are computed according to 

{0} = Z~ 1 d9i O[0] exp(-^W). (20) 

i 

0[6] denotes the dependence of the physical observable O on the configuration {0i}, the partition function 
Z is given by 

Z = J Y[d9i exp(-^W), (21) 

i 

where [3 = 1/ksT. The expectation values (20) are computed by means of the Monte-Carlo method using 
Wolffs 1-cluster algorithm [23]. This algorithm successfully tackles the problem of critical slowing down 
[24]. We computed the helicity modulus on lattices of different sizes L 2 x H, where H = 3,4,6,8,10 
and L = 40,60, 100 for each thickness. For some thicknesses we used L = 25,50, 100. Periodic boundary 
conditions were applied in all directions. We carried out of the order of 10,000 thermalization steps and 
of the order of 500,000 measurements. The calculations were performed on a heterogenes environment of 
workstations which include Sun, IBM RS/6000 and DEC alpha AXP workstations and on the Cray-YMP. 
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4 The two-dimensional region. 



Here we consider the temperature range T~ < T < T+ , where the 3D bulk correlation length exceeds the 
thickness of the film. This temperature range contains both the if-dependent 2D critical temperature T 2D 
as well as the 3D bulk critical temperature T\. For fixed H and for temperatures close enough to T 2D 
the system behaves as a 2D x — y model, thus, our method of analysis applied to the 2D x — y model [25] 
can be used here as well. In [25] we investigated the dimensionless quantitiy T/J, in the 3D system the 
helicity modulus aquires the additional unit of a length -1 , thus the proper quantity to consider is the ratio 
T(L,H,T)H/J. 

In the following sections we always leave H fixed. 

4.1 Finite-size scaling with respect to L above T^ D (H). 




Figure 2: The helicity modulus T(L, H, T) as a function of T for various lattices L 2 x 6. 

In Fig. 2 we show the data for T(L, H,T)a/J for films of fixed thickness H = 6 for various sizes L, a 
denotes the lattice spacing. As in [25] we can obtain a function T = F(L) via the beta function such, that 
T(L' , H , F(L')) = T(L, H , F(L)) = T p , where T p is a physical value of the helicity modulus. We define the 
beta function as: 

Since we expect our model to behave effectively two-dimensionally we use the ansatz 

P(T>T C (H)) = c(H)(T-T? D (H)) 1+ \ 

I3(T<T C (H)) = 0, (23) 
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Figure 3: The beta function obtained from the 100 2 x 6 and 50 2 x 6 lattices. 
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Figure 4: The helicity modulus T as a function of z for the 100 2 x 6 and 50 2 x 6 lattices. 



which is suggested by the Kosterlitz-Thouless theory. Inserting the ansatz (23) into (22) and integrating 
yields: 



InL 



B 

(T-T C (H)Y 



(24) 



where B = l/(i>c) and z is a constant of integration, which depends on the value of T used to define the 
scaling transformation. As in the pure 2D case we expect all values of T(L, H , T)/ J for fixed H to collapse 
on the same universal curve if T/J is considered as a function of z. As was explained in [25] this also means, 
that the correlation length £ grows according to £(T) oc exp B / (T — T c (H)y because T/J for fixed H should 
be a function of L/£(T) only, and thus, also a function of z = hi(L/£(T)). The beta function can be found 
numerically using methods described in [25]- [28]. 

Since the critical temperature T 2D , which enters in the expression for the beta function, is very sensitive 
to the planar extension L of the lattice [25] in contrast to the other parameters c and v, we only use the two 
largest lattices to numerically derive the beta function for each different thickness H . This way we determine 
a lower bound for the critical temperature T 2D (H). Having extracted data points for the beta funtion we 
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Table 1: Fitted values of the parameters c(H) and T c (H) of the beta function (23) for the used lattice 
pairs for different thicknesses H,v = 0.5. \ 2 an d the goodness of the fit Q are also given. 



fit them to the functional form (23) setting v = 0.5. The results of our fits are given in Table 1. Fig. 3 shows 
the beta function and Fig. 4 demonstrates that the values for the quantity Ta/J for the 50 2 x 6 and 100 2 x 6 
lattices collapse on one curve. 



4.2 Finite-size scaling with respect to L below T^ D (H). 

Since our calculation was performed on finite lattices of size L 2 x H we wish to reach the L — ► oo limit and to 
obtain T 2D (H) and b(H) (see Eq. (7)). In order to do that we need to know the leading finite- L corrections 
to the dimensionless ratio K = T/(TH). As the system behaves effectively two-dimensionally we can apply 
the formulae derived in [25] for the ratio K. This ratio satisfies the renormalization group equations (5) 
and (6), thus solving these equations for a finite scale / = ln_L in the limit L — ► oo and close to the critical 
temperature T 2D (H) yields [25]: 

K(L^oo,T<T c (H)) = Koo (T) (l + ^7^^% ) , (25) 

K(L,T cm = f(l-^), (26) 

where Koo(T) = K(L — ► oo,T) and 2 and 5' are constants of integration. 

At a fixed thickness H below T? D (H) we can use the expression (25) to extrapolate our values for K(L, T) 
at finite L to the values Koo(T) at infinite L. In order to do that we fit our calculated values for K(L,T) 
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Table 2: Fitted values of the parameters 2 and K^^H, T) of the expression (25) for various temperatures T 
and different thicknesses H . 



to the functional form (25). Table 2 contains our fitting results and Fig. 5 shows a typical fit. Note that the 
if-dependence of the ratio K is contained in the if-dependence of the critical temperature only. 
The extrapolated values for Kco{T) at a fixed H should behave as [18]: 
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(27) 



where we introduced the expected if-dependence of the parameters b and T^ D . At a fixed H both parameters 
can be determined by a fit of our results for Kco{T) to the functional form (27). The results of our fits are 
presented in Table 3. In Fig. 6 the fit to the data for A'oo(T, H) at H = 6 is shown. 
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Table 3: Fitted values of the parameters b(H) and T^ D (H) of the expression (27) for different thicknesses 
H . x 2 an d the goodness of the fit Q are also given. 
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Note that the if-dependence of the ratio K is contained in the if-dependence of the critical temperature 
T^ D and the constant b. These two parameters are the only parameters which are free to depend on the film 
thickness if we request the form of the equations (25), (26) and (27) to be applicable to our case of a film 
geometry. 

We would like to mention one difficulty in extracting the infinite L results (T) and finding the 
parameters T^ D (H) and b(H). In our study of the pure 2D x — y model [25] we found that formulae (25) 
and (27) are valid in a narrow interval below T^ D , in particular the expression (27) is a good approximation 
for temperatures in the region ~ 0.9T? D <T< T? D . For films of finite thickness H we have the additional 
problem of a crossover between 3D and 2D behavior, which makes the region of validity of expression (27) 
even smaller. Thus, one has to be careful in chosing the correct temperature range for the extrapolation 
procedure. 

5 Scaling of T/(TH) with respect to H. 
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Figure 7: T/(T(H,T)H) as a function of tH 1 !". v = 0.6705. Only even number thicknesses are displayed. 

This section is concerned with the scaling relation (10). In order to check the validity of the scaling form 
(10) we plot first T/(TH) versus tH 1 !" using the experimental value of Goldner and Ahlers v = .6705 [29] 
and the result of reference [13] T\ = 2.2017. This is done in Fig. 7. In the plot we only use temperatures 
below T* D {H). If the scaling beh avior (10) is valid, all our data points should lie on one universal curve. 
This is not the case for the films of thickness H = 3,4. For the other films scaling seems valid in the interval 
1.0 < tH 1 !" < 2.5. All the represented data have lost their Independence within errorbars. Scaling is 
confirmed by Fig. 8 which shows the same plot as Fig. 7, but only the data for the thicknesses H = 6,8, 10 
are displayed. The data points collapse onto one universal curve in the specified interval. The films with 
thickness H = 3,4 are too thin in order to exhibit the scaling behavior (10). 
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Figure 8: T/(T(H,T)H) as a function of tH 1 /". v = 0.6705. The solid curve is expression (17) with the 
parameters (28). 

We would like to compare our findings to the experimental results for the superfluid density of very thick 
films by Rhee, Gasparini and Bishop (RGB) [3]. In our language they plot TH versus tH 1 !" and find that 
their data do not collapse for the expected value of v. RGB demonstrated the lack of scaling of their data 
by collapsing them on one universal curve using a different value of v (v = 1.14(02)). Since we confirm 
scaling for the ratio K(T, H) in our simulation the lack of scaling of the experimental data is not due to a 
breakdown of the phenomenological scaling theory. Note that our results have been obtained with periodic 
boundary conditions. We are in the process of repeating the calculation using Dirichlet boundary conditions 
which we believe represent the experimental situation more closely. 

Since the scaling function <&(x) is known for temperatures which correspond to values of x close to x c we 
can actually find the constants A and x c in expression (17). For this purpose we plot K^^H, T) given in 
Table 2 versus tH 1 !" with the values for T\ and v as above for the three thickest films and fit the resulting 
data points to the form (17). We obtain 

A = 0.5928 ± 0.0050, 

x c = 0.9964 ± 0.0023. (28) 

The fitted curve is the solid line in Fig. 8. The universal function (17) with the parameters (28) describes 
the collapsed data rather well in the interval 1.0 < tH 1 !" < 1.3. Since T/T(T) oc £t{T) we would like to 
point out that $(1.3) ~ 1, i.e. £t ~ H . This agrees with the general picture that the 2D behavior sets in 
when the 3D bulk correlation length exceeds the film thickness. 
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6 independence of fitting parameters. 



In this section we would like to examine the if-dependence of T c and b. 
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Table 4: Fitted values of the parameters x c T\,l/v, and T\ of the expression (12) for various numbers of 
datapoints included in the fit. \ 2 an d the goodness of the fit Q are also given. 
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Figure 9: T^ D (H) as a function of H . The solid curve is the fit to (12), the dashed line represents (29). 

We fit our results for T? D (H) (see Table 3) to the expression (12) as expected. The results of the fit to the 
form (12) for all critical temperatures and the critical temperatures corresponding to the four thickest films 
are collected in Table 4, respectively. The fit including the critical temperatures of the four thickest films is 
shown in Fig. 9. The bulk critical temperature T\ obtained from the fit (which uses five data points) agrees 
within error bars with T\ = 2.2017 ± 0.0005 obtained by Janke [13] and Hasenbusch and Meyer [12]. The 
critical exponent turns out to be v = 0.6892 ± 0.0057 which is a little larger than the expected value. The 
same fit to only four critical temperatures yields v = 0.657 ± 0.012 and T\ = 2.1947 ± 0.0025. The critical 
exponent agrees with the expected value within errorbars, but T\ is somewhat smaller than Janke's value. 
The critical exponent v was estimated experimentally as well. In experiments on bulk helium Goldner and 
Ahlers [29] and Singsaas and Ahlers [30] deduced v = 0.6705 ± 0.0006 and v = 0.6717 ± 0.0004, respectively. 

We can obtain the same curves as described above using the parameters (28) entering the scaling expres- 
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Figure 10: ln&(if) as a function of In if. The solid curve represents (30). 



sions (12) and (16). The dashed line in Fig. 9 represents the function 

0.9964 

^yl./0.6705 



Tf u {E) = 2.2017 1.0 



(29) 



This curve is in fairly good agreement with the critical temperatures for the three thickest films, which we 
should expect because of scaling for these films. In Fig. 10 we plot b versus H and the function 



b(H) = 0.5928F 5/0 6705 . 



(30) 



Here the agreement for the three thickest films is not as good as for the case of the critical temperatures 
T^ D (H). However we attribute the scattering of the data points around the curve (30) to the procedure 
we applied in order to extract b(H) and T^ D (H). A possible cause of the poor agreement may be that 
more data points in the 2D critical region are required. Furthermore we have always assumed in the fitting 
procedure that the square root form in the expressions (7) and (17) is correct. A general exponent /j, leads 
to the following scaling form for b: 



b(H) = AH"!". 



(31) 



A slight deviation of /j, from 0.5 could also account for a better agreement of the curve (30) with the results 
for b(H). 

In summary we can state that our results for the critical temperatures T^ D (H) and for the parameter 
b(H) are in agreement with the expected ff-dependence, which is deduced from general scaling ideas. 
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7 Summary 

We have investigated the finite-size scaling properties of the superfluid density of a superfluid with respect 
to the film thickness. This was done by means of a Monte-Carlo simulation of the x — y model in a 
L x L x H geometry with periodic boundary conditions in all directions. We extrapolated the values of the 
superfluid density to the L — ► oo limit in the critical region where the model is effectively two-dimensional 
using the Kosterlitz-Thouless-Nelson renormalization group equations (5) and (6). The test of the scaling 
expression (10) revealed that scaling for the quantity T/(T(T, H)H) is fulfilled, i.e. the numerical results 
for this ratio collapse on one universal curve for sizes of film thickness H = 6,8, 10 used in our Monte- 
Carlo simulation. Furthermore we derived an analytic expression for the universal curve, which is valid for 
temperatures close to the critical temperatures of each film. Using the expression (7) we were able to extract 
the critical temperatures T^ D (H) and the parameters b(H) entering Eq. (27). The if -dependence of these 
two parameters agrees with the expected behavior (12) and (16) deduced from general scaling arguments. 
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